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We present a stochastic description of a model of N mutually interacting active particles in the 
presence of external helds and characterize its steady state behavior in the absence of currents. To 
reproduce the effects of the experimentally observed persistence of the trajectories of the active 
particles we consider a Gaussian forcing having a non vanishing correlation time r, whose finiteness 
is a measure of the activity of the system. With these ingredients we show that it is possible to 
develop a statistical mechanical approach similar to the one employed in the study of equilibrium 
liquids and to obtain the explicit form of the many-particle distribution function by means of the 
multidimensional unified colored noise approximation. Such a distribution plays a role analogous 
to the Gibbs distribution in equilibrium statistical mechanics and provides a complete information 
about the microscopic state of the system. From here we develop a method to determine the one 
and two-particle distribution functions in the spirit of the Born-Green-Yvon (BGY) equations of 
equilibrium statistical mechanics. The resulting equations which contain extra-correlations induced 
by the activity allow to determine the stationary density profiles in the presence of external helds, the 
pair correlations and the pressure of active huids. In the low density regime we obtain the effective 
pair potential ())(r) acting between two isolated particles separated by a distance, r, showing the 
existence an effective attraction between them induced by activity. Based on these results, in the 
second half of the paper we propose a mean held theory as an approach simpler than the BGY 
hierarchy and use it to derive a van der Waals expression of the equation of state. 


I. INTRODUCTION 

The collective behavior of microscopic living organisms and active particles capable of transforming chemical energy 
into motion has recently attracted the attention of the soft-condensed matter community as they present striking 
analogies, but also intriguing differences with respect to colloidal and molecular fluidJi®. Can we understand their 
rich phenomenology by applying experimental, theoretical and numerical techniques which proved to be successful in 
condensed matter physics? Recent investigations based on physical models akin to those widely employed in statistical 
mechanics, thermodynamics and rheology seem to give support to this hypothesis. 

Among the most studied systems we mention swimming bacteria, colloidal particles immersed in bacterial suspen¬ 
sions and self-propelled Janus particles. Swimming bacteria can be schematized as objects moving with speed vq along 
straight paths of average duration t, the so called persistence time, after which a random reorientation (a tumble) 
takes place. After many such events, i.e. on a time scale larger than r, the particle displacement can be assimilated 
to a random walk and diffusive behavior emerges with a diffusion constant, Da = v^t/ 2d, where d is the space di¬ 
mensionality. The resulting motion displays two peculiar features: a) the trajectories display an anomalously long 
persistence not observed in Brownian motion, i.e. their direction and the velocity remain constant for a time lapse 
much longer than those corresponding to colloidal particles, and the similarity with ordinary diffusion appears only 
on a longer time-scale; b) there is a spontaneous tendency of the particles to aggregate into clusters notwithstanding 
there is no evidence of direct attractive forces, while on the contrary short repulsive inter-particle forces are at work. 
This is in a nutshell the idea of the successful run-and-tumble (RnT) model that captures many aspects experimentally 
observecPlt^. Experiments conducted employing swimming bacterial suspensions have shown that their diffusivity can 
be hundred times larger than the one arising from the thermal agitatiorP®^. To give some numbers the diffusivity 
of Escherichia coli bacteria is D ^ 100 /s, whereas passive bacteria diffuse with D « 0.3 pLm?/s. Active Janus 

particles instead have ZJ in a range 4 — 25 jjLm?/s, finally colloidal particles immersed in a bacterial suspension display 
D « 0.5 — 0.1 firri^/s. In all cases, the contribution due to the diffusion due to the thermal agitation of the molecules 
of the solvent is more than ten times smaller than the one due to the activitjE^. 

The work of several groups has led to the formulation of theoretical models and to the application of methods of 
statistical mechanic^2®3 whereas Brady and TakatorP^ have put forward an approach based on thermodynamics 
and microrheology and introduced the seminal idea of swim temperature Ts and swim pressure T^p, p being the 
particle number density. In addition to the swim pressure, the total pressure contains other contributions stemming 
not only from the particle excluded volume but also from the effective attraction experienced by the active particles. 
Such an effective attraction is a peculiar aspect of active matter: motile particles, characterized by isotropic repulsive 
interactions due to the persistence of their motions become slower at high density and tend to form clusters and/or 
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pile-up in regions adjacent repulsive substrates. The attractive forces may eventually lead to a scenario where van der 
Waals loops appear as soon as the off-equilibrium control parameter, the persistence time, is above a certain critical 
value. 

The idea of the present paper is to use a microscopic model capable of reproducing the basic features of the 
RnT model and of the active Brownian particles. In our model the microscopic state is specified only by the particles 
positions, while the effect of the angular dynamics is encapsulated in a Gaussian colored noise having a finite relaxation 
time, reflecting the persistence of the directions. The correspondence between the models can be established by 
noting that, in the free particle case, in all models each velocity component has an exponential time-autocorrelation 
function characterized by a correlation time Despite the differences in the driving stochastic forces, as shown 

in Ref.l^, Gaussian-colored noise driven particles behave strikingly similarly to RnT particles when in presence of 
steeply repulsive interactions. Moreover the Gaussian colored noise model is the only one for which an approximate 
stationary distribution for multiple particles is knowrl^. This motivates us further in studying the Gaussian colored 
noise model in more detail. 

Somehow surprisingly the present work shows that the approximate stationary configurational distribution of a 
system of interacting particles undergoing over-damped moti on un der the action of Gaussian colored noise bears 
a strong resemblance with the equilibrium Gibbs distributiorP^, with the due differences: a new temperature, 
named swim temperature, takes the role of the ordinary temperature, the probability of configurations depends on a 
complicated function of lA, the potential energy. Such a result is independent on the fine details of the interactions, 
but certainly depends on the persistence of the random fluctuations. The most direct consequence of the form of 
this steady state distribution is the appearance of an effective attraction between two or more active particles in the 
absence of any attractive direct interactiorP^. The reduction of the mobility of the particles due to the presence of 
other particles in their vicinity may eventually lead to a phase se para tion between a high and a rarefied density phase, 
a phenomenon named motility induced phase separation (MlPSjP^. 

This paper is organized as follows: in Sec. |TT]we present the coarse grained stochastic model describing an assembly 
of active particles, consisting of a set of coupled Langevin equations for the coordinates of the particles subject to 
colored Gaussian noise. We switch from the Langevin description to the corresponding Fokker-Planck equation for 
the joint probab ility distribution of N particles and within the multidimensional unified colored noise approximation 
(MUGNAjpil^, in the stationary case and under the condition of vanishing currents we obtain its exact form. The 
obtained distribution implies detailed balance and potential conditions as discussed in appendix G. In the case of 
non vanishing currents it is straightforward to write such a distribution for a single particle in one dimension in the 
presence of colored nois^^, whereas in higher dimensions and in interacting systems sucha a generalization cannot be 
obtained by the method of appendix A. Using the 7V-particle distribution in Sec. |III| we derive the first two members 
of a hierarchy of equations for the marginalized probability distributions of one and two particles that play the role 
of the Born-Green-Yvon (BGY) equations for active systems. We apply the first of these equations to the study of a 
density profile in the presence of an external potential, extend the treatment to the case of interacting active particles. 
In the low density limit we are able to write the exact form of such the pair correlation and define the effective pair 
potential between two isolated particles. In section |IV| we employ a variational approach, complementary to that 
based on the BGY equations which are exact but highly unpractical, and taking advantage of the explicit form of 
the MUGNA equation we construct a "free energy" functional whose minimum corresponds to the exact solution. 
The mean-field theory for interacting systems follows by searching the solution among the probability distributions 
which are a product of single particle probability distributions. The functional is finally used in section |V] to lay out a 
method capable to interpret the phase behavior of the model in terms of the relevant control parameters starting from 
a microscopic description. Finally, we present our conclusions in |VI[ We include with four appendices: in appendix [A| 
containing the calculation details leading to the MUGNA, in appendix]^ we discuss the approach to the solution and 
establish an H-theorem, in appendixwe verify the detailed balance condition, whereas in appendix]^ we discuss a 
technical aspect which allows to evaluate a key ingredient of our approach, the determinant of the friction matrix. 


II. THE MODEL AND ITS STATIONARY MANY-PARTICLE DISTRIBUTION FUNCTION 

In order to describe the properties of a suspension of active particles we consider a three dimensional container of 
volume V where an assembly of N interacting active spherical particles at positions r^, with i ranging from 1 to N, 
subject to external fields undergoes over-damped motion driven by random fluctuating forces of different origin and 
nature. In fact, each particle besides experiencing a white noise force, due to the thermal agitation of the solvent and 
characterized by a diffusivity Dt , is acted upon by a drag force proportional to its velocity, r^, and in addition is 
subject to a colored Gaussian noise, Vi, of zero mean, characteristic time r and diffusivity Da- This second type of 
noise is intended to mimic on time scales larger than r the behavior of self-propelled particles whose propulsion force 
randomizes with characteristic time r. Such a model involving only the positional degrees of freedom of the particles 
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has the advantage of allowing for analytical progress and for this reason can be used instead of more microscopic 
models where the rotational dynamics is fully accounted for. We consider the following set of equatio ns of m otion for 
the positions of the active Brownian spheres which has been treated in the literature by some author^l ^^ l ^^ l : 


dt) = ... ,rN) + + Vi(t) 
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are coupled to changes of the velocities v described by: 


1 

v^{t) = - Vi{t) + 

T T 


( 1 ) 


( 2 ) 


The force Fi = —ViU acting on the i-th particle is conservative and associated to the potential U{yi,. .. ,rjv), 7 is 
the drag coefficient, whereas the stochastic vectors and J 7 j(<) are Gaussian and Markovian processes distributed 
with zero mean and moments = 2SijS{t — t') and {'ni{t)riAt')) = 2SijS{t — t'). While Dt represents a 


translational diffusion coefficient, the coefficient Da 
Uhlenbeck process Viit) via 


due to the activity is related to the correlation of the Ornstein- 


{v,{t)vj{t')) = 2 — 6ij exp(-|f - t'l/r). 

T 

In spite of the fact that the magnitude of the velocity has not the fixed value uq as in the RnT, but fluctuates, a 
correspondence can be established between the present model and the RnT by requiring that the mean value (vf) = Vq, 
so that Da = VQT/2d. 

In the following we shall adopt a shorthand notation in which the set of vectors is represented by an array {cci} 
of dimension M = dN and the remaining terms are replaced by non bold letters. We assume that the force Fi may 
be due to action of external agents and mutual interactions between the particles. In order to proceed we consider 
the multidimensional version of the unified colored noise approximation (UCNAj^ which consists in eliminating 
adiabatically the fast degrees of freedom of the problem and, as shown in appendix]^ arrive at the following equation 
for the the particles coordinates: 






L 7 


( 3 ) 


where we introduced the non dimensional friction matrix Tik 


Bifc — S’, 


'ik 


T d^u 

7 dxidxk' 


( 4 ) 


Such a formula says that the effective dynamics of each particle depend on its distance relative to the other particles 
and on its absolute position if an external field is present, not only through the direct coupling Fk, but also through 
the motility matrix, ^T^^ which is the sum of a constant contribution due to the background fluid plus a space 
dependent term due to the interparticle forces mediated by the colored bath. 

Notice that such a structure of the friction matrix introduces velocity correlations among different velocity com¬ 
ponents of a given particle or between the velocities of different particles (this aspect will be discussed in detail in a 
forthcoming publication). For this reason, the present approach is not a mapping onto a passive equilibrium system, 
but rather represents a mapping onto a system with a self generated inhomogeneous friction. This leads to strong 
deviations from equilibrium such as the explicit dependence of the stationary state on the transport coefficient. 

In order to obtain meaningful results one must require that all eigenvalues of T are non negative. We cannot prove 
such a condition in general, however, it seems a reasonable assumption for repulsive pair potentials. On the contrary, 
it is easy to find examples where appropriately chosen external potentials determine negative eigenvalues and limit 
the validity of our formula. The structure of eq. ^ together with Q is interesting because it tells that the damping 
experienced by a particle is due to a standard drag force —yx plus a contribution stemming from interactions. Thus 
the effective friction increases wit h de nsity leading to lower mobility and to a tendency to cluster. This mechanism 
can eventually lead to the (MIPS)I^I^ and is an intrinsically non-equilibrium effect. 

Let us write the Fokker-Planck equation (FPE) for the iV-particle distribution function /at associated with the 
stochastic differential equation eq. ([^ under the Stratonovitch conventiorP^ 


dfNixi 


, Xat , t) 


d 


dt 


- = - Ji(xi,...,XAr;t) 


( 5 ) 
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where the l-th component of the probability current is: 


^ Err.'- (^a + A) E A]) 


( 6 ) 


Using the method illustrated in Appendix and requiring the vanishing of all components of the probability current 
vector Ji, without further approximations we obtain the following set of conditions for the existence of the steady 
state iV-particle distribution function Pjy: 




Pn E ^— 


T d^U \dU{xi,...,XN) 


7 dxndxk! dx 


= 0 . 


(7) 


We now define the effective configurational energy of the system T-L{xi,... ,Xn) related to the bare potential energy 
li{xi,.. .,Xn), by: 


N 

'H{xi,...,xn) =U{xi,...,xn) + ^E( 


dU{xi,.. .,xn)\'^ 


dxk 


) -(£>a + A)7ln|detrifc| 


( 8 ) 


It is easy to verify that the following iV-particles configurational probability distribution, obtained for the first time 
in ref.^, 


Pn{xi, ■ ■ ■,xn) = E expf- 
Zn V 


Hixi,.. .,xn) \ 

{Da + Dt)-f ) 


(9) 


is a solution of eq. 0 , where Zjv is a normalization constant, the analogue of the canonical partition function and 
enforcing the condition f ... f Pn{xi, ..., xn) dxi... cLxn = 1- Notice that the mobility enters the stationary 
distribution, at variance with equilibrium systems where the kinetic coefficients never influence the form of the 
probability distribution function. In the limit r —>■ 0, Zjq reduces to the equilibrium configurational partition function 
for a system characterized by energy U. Formula 0 is a generalization to spaces of arbitrary dimensions of a stationary 
distribution obtained by Hanggi and Jung for a single degree of freedom in one dimensioipSl. It displays an exponential 
dependence on a function constructed with the potential lA and its derivatives. These derivatives must be non singular 
and also satisfy the condition that the determinant is non negative. In the case of a single particle the result can 
also be derived by a multiple-time scale method often employed to reduce the phase-space Kramers equation to the 
configurational space Smoluchowski equatiorP^. Notice, that the explicit form of the distribution allows to identify 
an effective temperature of the system with the quantity 


T, = (A + A)7, (10) 

which shows the typical Einstein fluctuation-dissipation relation (with = 1) between temperature, diffusivity 
and and drag coefficient. Let us remark that the validity of formula ([^ is limited to the regime Da > Dt, which 
also corresponds to the actual experiments with bacterial baths, whereas the limit Da —t 0 with r and Dt finite, 
neither corresponds to the real situation but also leads to meaningless theoretical predictions since it violates the 
hypothesis under which the MUCNA is obtained. To conclude this section, we identify the temperature with Tg with 
swim temperatur^^ and introduce as a measure of the distance from thermodynamic equilibrium the Peclet number, 
Pe = y/D^/ a, which is the ratio between the mean square diffusive displacement due to the active bath in a time 
interval r and the typical size of the active particles, say a. 


III. BGY HIERARCHY AND FLUID STRUCTURE 

As it stands formula ^ is exact, but contains too many details to be really useful, however, borrowing methods of 
equilibrium statistical mechanics we can trace out degrees of freedom and arrive at formulas involving the distribution 
functions of only few particles. To this purpose we shall use the stationarity condition Q to derive a set of equations 
equivalent to the BGY hierarchy for equilibrium correlations. The hierarchy becomes of practical utility if utilized in 
conjunction with a suitable truncation scheme in order to eliminate the dependence from the higher order correlations. 
Let us turn to standard vector notation where the indices a and j3 running from 1 to d identify the Cartesian 
components and latin subscripts the particles. The total potential is assumed to be the sum of the mutual pairwise 
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interactions w{r — r') between the particles and of the potential exerted by the external field u(r): U{yi^ ... = 

w{Yi,Yj) + u{Yi). The hierarchy follows from writing eq. 0 in the equivalent form: 


H • ■ • ’ rjv)-Piv(ri,..., yn)] = Pn^Vi, 

/3 ra k^l 


dw{Yi - Yk) 

dr„i 


( 11 ) 


We proceed to marginalize the N dimensional distribution function P/v introducing the reduced probability distri¬ 
bution functions of order n as Pj^\Yi,X 2 , ■ ■ ■, fn) = / dr„+i... drArPAr(ri, r 2 ,..., yn). By integrating eq. Q over 
(N — 2) coordinates we obtain an equation for pj^'^ (ri, r 2 ) in terms of higher order marginal distributions and choosing 
n = 2 we find: 


-T. 11 dr ,.. .*» EE ■ ■ -.r™)] = Pg'{r,.r,){^^ + 

+ 51 / drfcP®(ri,r 2 ,rfc) 


/3 n 

.(3). „ ^,dw{Yi-Yk) 


k>2 ‘ 


dr, 


( 12 ) 


Now, we notice that in the case of a large number of particles and in the limit of small r/y the matrix T^j*^ is nearly 
diagonal and can be approximated by 


n-l 

al 




k^l 


where Ua /3 = and Wa /3 = gr ^drl • Substituting this approximation in eq. (12) we find: 




dr pi L 


N 


{Yi,Y2)Sap - ^(pt^lri,r2)u„/3(ri) -h Pt^^(ri,r2)Wa/3(ri - r2) + J '^dYkPj:^\Yi,Y2,Yk)Wap{Yi - rfe)^ 


= -p. 


( 2 ) 


N 


(l'l,l'2)( 


5u(rai) 5u;(ri-r2) 


Oral 


+ 


dr, 


al 


- E / (ri,i- 2 ,rfc) 


k >2 ‘ 


,dw {Yi - Yk) 
Oral 


(13) 


which represents the BGY equation for the pair probability distribution P 2 . By integrating also over the coordinate 2 
and switching to the n-th order density distributions: (ri, r 2 ,..., r„) = r 2 ,..., r„) we obtain the 

first BGY-like equation: 




d 


dr pi 1 


Sa/ 3 P^^\ri) - -p^'^\Yi)UapiYi) - - J cIy 2 {y i, Y 2) Wap{Yi - Y2) 


drai J dr^i 


(14) 


that in the limit of r —0 is just the BGY equation for the single-particle distribution function. By performing an 
analogous substitution in (131 we obtain the BGY equation for the pair correlation function including the corrections 
of order r /7 stemming from the activity. 

The r.h.s. of equation (141 contains the coupling to the external field and the so-called direct interaction among the 
particles, whereas the l.h.s. besides the ideal gas term contains a term proportional to the ac tivity parameter that we 
name indirect interaction term, following the nomenclature introduced by Solon et al.l^^l^, although our expression 
does not coincide with theirs because the present model does not depend on angular degrees of freedom. Notice that 
here the indirect interaction term stems from the expansion of the determinant to first order in the parameter r/y. 
On the other hand, if one employed an higher order expansion in this parameter terms involving terms up to Wbody 
correlations would appear. 


A. Non interacting active particles under in-homogeneous conditions 


In the case of vanishing "inter-molecular" forces equation (141 gives access to the single particle distribution of 


Brownian active particles near a wall. Let us assume w{r) = 0 and u(r) = u{x) to represent a generic (flat) wall 
confining potential, twice differentiable and with the properties that lim 2 :_,._oo u{x) = 00 and lima;_>oo u{x) = 0 . Since 
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the particles are non interacting it is not necessary to expand the friction matrix in powers of activity parameter as 
in eq. (141 and one can use its exact expression: 

^ Sin- 




Alternatively, one can regard such a formula as a resummation of the generalized binomial series (1 — y + ...), 

whose first term, y = -Uxx{x) appears in the l.h.s. of ([l4|. We write: 


^^dxil + lUxlix)) ~ 


and find the profile: 


p^^\x) = poexp - 


+ i:iiux{x)f 


1 -f Uxx {x^ 
7 


(15) 


(16) 


Using the mechanical definition of pressurd^SEl this simple example provides an exact expression for the pressure 
exerted by non interacting active particles: in fact, by integrating both sides of (151 with respect to x from —oo to 
oo and recalling that the right hand side is just the force per unit area exerted by the wall, located at a: « 0 , on the 
fluid we obtain the pressure as: 


Ps = Tsp{x oo) 


(17) 


having assumed that the negative region is not accessible to the particles. Notice that Ps = pTg does not depend 
on the particular form of the wall potentiaP^. Such a pressure is the so called swim pressure that is the sum of the 
active and passive ideal pressures. As remarked by Brady ps may depend on the size of the particles only through 
the hydrodynamic drag factor 7 and in the present case contains a thermal contribution associated with Dt and an 
athermal part associated with the active dynamic^^. 

Repulsive barriers with positive curvature {uxx > 0) induce a local accumulation of particles and lower their mobility 
with respect to their bulk value. One can speculate that a similar phenomenon occurs spontaneously in an interacting 
system where denser and less motile clusters of active particles "attract" fast moving particles from rarefied regions: 
the flux would be sustained by the difference between the pressures of the two regions. We shall consider the role of 
interactions among the particles in the section below. 

Interestingly, we can rewrite (151 as the local balance equation between a local pressure term Ps{x) = Ts{x) p^^\x) 
and a force term with the local swim temperature Ts{x) defined as 


fs{x) = 


T, 


1 “t” Uxx 


{x)' 


(18) 


In Fig. [T] we display both the density and temperature profiles next to a repulsive wall for three different values of 
the persistence time. Formula (151 applies to rather general potentials under the condition that Uxx{x) > —q/r. 

It is of interest to apply (15) to the special case of sedimentation of active colloidal suspensions, where u{x) = 
mgx. One immediately finds the barometric law p(x) = po exp{—mgx/Ts), predicted by Tailleur and Cates in the 
small sedimentation rate limitP^ and confirmed experimentally by Palacci et al.^ who performed the so-called Perrin 
sedimentation experiment using active particles. This scenario was also confirmed for swimming bacteria under 
centrifugatioiPS. Notice that in this case the linear potential only appears in the exponent as if the system were at 
equilibrium, but with an effective temperature Tg higher then the ambient temperature Tg = Dt^ according to the 
formula Ts/Tq = 1 + Da/Dt- Formula (151 is more general and encapsulates the idea that repulsive interactions render 
the diffusion of the particles less efficient. 


B. Interacting active particles and non ideal pressure via BGY equation 


To extend eq.(15l to the case of interacting active particles we apply an external potential varying only in the 
direction x and factorize the multiparticle distributions as p^"'^(ri,..., r„) Ki p^^\ri)... p^^'> 
field-like approximation and recast (141 as: 


i) invoking a mean 


_T — f P^^\xi) 

®da;i yi + /:^Uxi,xi{xi) + ^ J dT2p^^'> {r2)wxi,xi{xi,r2)) 

p^'^\xi)^u{xi) + p^^\xi) [ dr2p^^\x2)-^w{ri,r2)) 
ciXi J UX\ 


( 19 ) 
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Density profile 



Figure 1: Three density profiles in the presence of a repulsive soft wall of the form u{x) = corresponding to three 

different values of the persistence time t as indicated in the box. In the inset we report the effective mobility profiles for the 
same values of r. 


where we have dropped off-diagonal terms Wap in the denominator by the symmetry of the planar problem. 

Formula (19) is valid up to linear order in r/y and in the l.h.s. contains the gradient of the density multiplied by 
a space dependent mobility, represented by the denominator. It generalizes to the interacting case the denominator 
featuring in formula (15) by the presence of inter-particle contributions. When Wxx(j) is positive it gives rise to a 
slowing down of particle motion due to the surrounding particles in interacting active systems and determines a density 
dependence of the mobility. In the r.h.s. the first two terms represent the external force term and the contribution 
to the pressure gradient due to direct "molecular" interactions, respectively. Notice again that the denominator 
represents the combined action of the "molecular" forces and the dynamics, and can be identified with the indirect 
interaction of Solon et al.. Eq. (19 1 in the limit t —>■ 0 reproduces the expression of the hydrostatic equilibrium of a 
standard molecular fluid in a local density approximation. There is an interesting analogy with granular gases subject 
to homogeneous shaking where the kinetic temperature can be a function of position and regions of higher density 
correspond to lower temperatures because the higher rate of inelastic collisions causes a higher energy dissipation 

ratepiU. 


C. Effective pair potential 


Equation (13) determines the pair correlation function of the model: /9*^^^(ri,r2) = p^g{\ri — r 2 |) , but it suffers of 


the usual problem of statistical mechanics of liquids since it requires the knowledge of higher order correlations even 
at linear order in the expansion in the activity parameter. One could proceed further by assigning a prescription to 
determine the three particle correlation function in terms of g{r) as done in the literature, but we shall not follow 
this approach and only consider the low density limit where it is possible to use it to derive a simple expression for 
g(r) and define an effective interaction. In fact, the pair distribution function for a two particle system obtained from 
^ reads: 


g{r) = exp - 


w{r) + ^[K(r)]2 - T, ln[(l + 2^w"{r)){l + 

T. 


( 20 ) 























where the apostrophes mean derivative with respect to the separation r. The effective pair potential is defined as 
= —Tg ln 5 (r). Let us remark that this method of introducing the effective potential does not account for the three 
body terms which instead would emerge from the solution of the BGY-like equation. The lack of such contributions 
affects the calculation of the virial terms beyond the second in the pressure equation of state discussed in section |V] 


IV. MEAN FIELD PSEUDO-HELMHOLTZ FUNCTIONAL 


As we have seen above, the BGY hierarchy is instructive, but unpractical because it requires a truncation scheme 
(for instance the Kirkwood superposition approximation at the level of two-particle distribution function to eliminate 
the three body distribution) and any progress can be obtained only numerically. The approximations involved are 
often difficult to assess, and for the sake of simplicity in the present paper we limit ourselves to a simpler mean field 
approach. A method often employed in equilibrium statistical mechanics to construct a mean field theory is the so 
called variational method based on the Gibbs-Bogoliubov inequalitjJ^. At equilibrium one can prove that there exists 
a Helmholtz free energy functional, of the probability density distribution, /, in configuration space such that 

it attains its minimum value, when the generic distribution /, selected among those which are normalized and non 
negative, corresponds to the equilibrium distribution function /®'^. Such a method can be generalized to the present 
non equilibrium case to develop a mean field theory. We start from the probability distribution (|^ and introduce the 
"effective Helmholtz free energy" functional as 

= TrP)([“' [n -k Tg In (21) 

where Tr = ff dri ,... ,drjv and is an arbitrary 7V-particle normalized probability distribution, thus integrable 

and non negative. Define 


Zpf = Tr exp 


'H(ri,... ,rAr) 


T, 


The stationary probability distribution which minimizes P is: 

Hixi,.. .,xn) 


In = exp 


Tg 




( 22 ) 


(23) 


and is the analogue of the equilibrium density distribution. In fact, the effective free energy: 

MfN] = -Ts\nZM] ( 24 ) 

is a lower bound for all distributions such that > Po[/^] for yk (see ref.l^. 

In other words, the functional evaluated with any approximate distribution has a "free energy" higher than the one 
corresponding to the exact distribution. Since T is minimal at the steady state it is then reasonable to assume that it 
can play the role of an Helmholtz free energy in an equilibrium system. Using such an analogy we shall construct an 
explicit (but approximate) mean field expression of F in terms of the one body distribution function Pi or equivalently 
Let us remark that our present results do not allow to establish the central achievement of density functional 
theory (DFT), namely the fact that the 7V-particle distribution function /(y is a unique functional of the one body 
density distribution^. This in turn means that for any fluid exposed to an arbitrary external potential u(r), the 
intrinsic free energy is a unique functional of equilibrium single-particle density . 

Nevertheless, it is possible to prove another remarkable property of P, an H-theorem, stating that the time derivative 
of P evaluated on a time dependent solution of the FPE is non positive, meaning that the functional always decreases 
during the dynamics and tends to a unique solution. The proof of this property is a simple application of a general 
result concerning the Fokker-Planck equation and can be found in the book by Risked^ and in appendix]^ we provide 
the necessary elements to obtain it. 

In order to construct the mean field approximation one assumes the following factorization of the Y-particles 
probability distribution, P^®“*(ri,... = n()(^j^Pi(rfc), and minimizes the functional 

p]^ptr^al^ = TrPj:^^^\n + Tg InP^™') ( 25 ) 


since the minimum of such a functional gives the closest value to the true functional PI/n] within those having the 
product form as above: 


= iVT, / drPi(r)lnPi(r)-krrP^™''H 


(26) 
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having used the fact that all Pi are normalized to 1. Switching to the distribution we find: 

N 

=tJ drp(i)(r)(lnp(i)(r) - 1) + ^(H„) (27) 

1 

where we considered the following decomposition of "H = T-Ln into one body, two body, three body up to N- 
body interactions and their averages (’H„) = J ■ ■ ■ J dri... (ri)... (r„)'H„(r,..., r„). The minimum of the 

functional P is obtained by differentiating w.r.t. to the p^^'>{r). 


A. Application to Planar and Spherical interfaces of non interacting systems 


As an example we consider the "Helmholtz functional" for non interacting particles near a flat wall, for which only 
Hi contributes: 


(X f dxp^^\x)(Ts{\n - — -l)-Ts ln[l + -u^x^x)] + u{x) + ^(ux(x) f) 

J ^ Po 7 27 / 


(28) 


minimizing w.r.t. p^^^ we obtain the same result as (161 which is exact. Similarly, the extension to three dimensional 
spherical walls is 


= [ drp^^\r)(Tsiln -—— - 1) - Tsln[(l + -ii"(r))(l + + u{r) + :^(u'(r))^) 

J \ Po 7 IX 2.'^ J 


(29) 


Notice the difference of order r/y, namely a factor 2, between the g(r) found by considering two active particles and 
the single-particle density profile when a spherical wall particle is pinned at the origin. Such a non-equivalence is 
fingerprint of the off-equilibrium nature of the system and disappears as r —> 0. 


1. Mobility and pressure tensor in spherical symmetry 


One can derive an interesting formula starting from (|29|). We first extremize the functional with respect to p^^\r) 

it bv nPl 


(i.e. SP/6p^^\r) = 0 ) to obtain the profile, hence differentiate the result with respect to r and multiply it by p' 
and arrive at the following formula: 



= -p(i)(rX„i(r). 


(30) 


On the other hand in a system of spherical symmetry the pressure tensor, pap(r), at point r must be of the form: 
Papir) = farppN{r) + (Sap — fafp)pT{r), where we have separated the normal (N) and tangential (T) components 
and r denotes unit vector in the radial direction. In the absence of currents the mechanical balance condition dictates 

^PN{r) P -{pn{t) - pT{r)) = -p(^)(rX„i(r). 
ar r 


By comparing with eq. (30l we may identify the two components pjv(r) = Tsp(r)/(1 -|- -u'^^tix)) pt{x) = 


Tsp^'^\r)/{1 -I- It would be interesting to extend the above analysis to the case of interacting particles 


along the lines sketched in section III , but this task is left for a future paper. Eq. (301 provides an interesting 


generalization of the planar formula (151 to a spherical interface and shows that unlike equilibrium fluids the normal 


and tangential components of the pressure tensor of a gas of non interacting active particles are not equal in the 
proximity of a spherical wall, but they tend to the common value Pp^^'^ when r —> 0 . Notice that pt{x) exceeds 
Pn{x) for a repulsive wall-potential since u'^xtip) < 0. By pursuing the analogy with equilibrium statistical mechanics 
one could identify the integral of the second term in eq. ( |30[ ) with a mechanical surface tension, which by our argument 
would turn up to be negative in agreement with Bialke et al.l^. A question arises: what is the nature of pjy and pr? 
The two quantities are two components of the swim pressure tensor and are different from their common bulk value 
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p, because near a a spherical obstacle the motilities of the particles along the normal and in the tangential plane to 
the surface are also different. As far as the mobility tensor is concerned in spherical geometry we can decompose it as 




1 


1 






Trir) 


{Sai3 - Tafp) 


(31) 


with rAr(r) = (1 + and Trir) = (1 + r)/r) which shows why the tangential motion is characterised 

by a higher mobility than the normal motion near a curved surface as reported by many authors on the basis of 
simulation results and phenomenological argument^^. Particles are free to slide along the directions tangential to the 
surface and this explains why px > Pn ■ 

Recently, Smallenburg and Loweiii^ have studied numerically non interacting active spheres in spherical geometry 
and found that for finite curvature radii active particles in contact with the inside of the boundary tend to spend 
larger time than those at the outside and differences between the inner and the outer density profile as a function of 
the normal distance to the wall. To see how the density profile around a spherical repulsive wall reduces to the profile 
near a planar wall we use the explicit form of the solutions and introduce the normal distance, z = r — Rq, to the 
spherical wall defined by the potential jir — i?o)" and write for z > Rq\ 


p(i)(z) = poexp 


1 

t; 



IH-yn(n + 

7 



27 (T^ 



2n+2' 


T Mo 1 
- n 

7 cr ito + 



(32) 


In the limit of z/i?o 0 the spherical profile must reduce to the planar profile, and indeed this is the situation as one 
can see by expanding the formula above in powers of z/Rq • For particles contained in a spherical cavity {z < Rq) 

the potential reads Mo(T"/(i?o ~ ’’)" and one must replace the last factor by ^1 + 7 ^■ "^he average 
local density on the concave side of a sphere is larger than the corre^onding density on the convex side as found 
numerically by Mallory et al.l^. This is illustrated in Fig. [^and Fig. 


Density profile 



Figure 2: Comparison between the outer and the inner local density profiles for three choices of r/7 in the presence of a 
repulsive spherical soft wall at Ro = 10. Particles are non interacting. 
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Pressure profiles 



Figure 3: Comparison between the outer and the inner local pressures prohles for two choices of r/7, 0.1 and 0.3 in the presence 
of a repulsive spherical soft wall at Ro = 10. The tangential component displays a peak at a distance ~ a from the wall, whereas 
the normal component decreases monotonically as the wall is approached. Notice that the tangential components are always 
larger than the corresponding normal components. 


B. Can we write a density functional for interacting active particles? 


As discussed above, since in the case of non interacting systems it is straightforward to construct an appropriate 
"Helmholtz" free energy functional such that it yields the same equation as the BGY method, one would like to 
determine the corresponding functional also in the case of interacting particles even within the simplest mean-field 
approximation. To this purpose we should construct a mean field functional T whose functional derivative with 
respect to gives a non equilibrium chemical potential whose gradient must vanish in the steady state 

giving a condition identical to the BGY equation (141. We have found that such a program can be carried out in the 


one dimensional case, in fact the following construction has the required properties 

J'[p(i)]=Ts J dxp'^^^x) {In - 1)-Ts^ J dxuxx{x)p^^\x) - jj dxdx'Wxx{x,x')p'^'^\x)p'^'^\x') 

+ j dxu{x)p^^\x)+ ^ JJ dxdx'w{x,x')p^^\x)p^^\x'). (33) 

Unfortunately, we are not able to write the corresponding functional for dimensions higher than one, the difficulty 
being the presence of the off-diagonal elements of the friction matrix feauturing in eq. (141 which render the equation 
non integrable. In other words, eq. (141 is only a mechanical balance condition and involves transport coefficients, 


such as 

coefficients never appear when thermodynamic equilibrium holds. 


Vij, thus revealing the true off-equilibrium nature of the active system, while in passive systems the transport 


V. VAN DER WAALS BULK EQUATION OF STATE FOR ACTIVE MATTER 


As we have just seen the construction of a free energy functional capable of describing the inhomogeneous properties 
of active particles remains an open problem, nevertheless one can focus on the less ambitious task of determining the 
pseudo-free energy of a bulk system characterized by an effective pair interaction (j){r), defined in Section III Notice 
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that in the homogeneous case the friction matrix reduces to a scalar quantity due to the higher symmetry. As a 
prerequisite for a successful mean-field approximation (f>(r) must be splitted into a short range repulsive part, (f>rep, 
and a weaker longer range attractive contribution, (patj and only the latter can be treated in a mean field fashion, since 
the effect of repulsion is a highly correlated phenomenon and not perturbative. Explicit expressions for the splitted 
potentials are given by (361 and (37). To capture the effect of repulsive forces a simple modification of the ideal gas 


entropic term, already introduced by van der Waals to account for the reduction of configurational entropy due to the 
finite volume of the particles, is sufficient. Thermodynamic perturbation theory would represent the natural choice to 
determine the total free energy. It assumes a reference hard-sphere system, whose Hamiltonian only depends on (f>rep, 
as the unperturbed system. In particular the reference system characterized by a temperature dependent diameter, 
d{T), allows to determine the free energy excess associated with the perturbing potential (fiat and and equation of 
state of the model. However, since perturbation theory still requires a considerable amount of computer calculations, 
here we make the simplest ansatz and write the following free energy functional 


^vdW 


r) 


In 




l-fp(i)(r)d3(r,) 


-1 


^ J J\r'-r\>d 


The first term is just the entropy of a fluid with the excluded volume correction, the second term stems from the 
activity and may lead to condensation phenomena for large values of r. In relation to the naive mean-field functional 
(331, the van der Waals entropic term already contains the repulsive part of the direct interaction, whereas the 
attractive term which vanishes in the r —>■ 0 limit takes into account the activity. In principle this form of free 
energy can be constructed by employing the Mayer cluster expansion to evaluate Zn, associated with the effective 
Hamiltonian 77, in the approximation of neglecting m-body interactions with m > 2 there contained, that is using the 
effective potential (f>{r). With these ingredients one can define the active component of the van der Waals attractive 
parameter: 


a(A) = _(i/2) [ d^r^atir) 

Jr>d{T,) 


(35) 


and the covolume 


b = 


2tt 


d^{T,) 

3 


and estimate J^vdw- Since the procedure adopted is completely analogous to the one employed in equilibrium fluids 
one can immediately derive a pressure equation by differentating !Fvdw with respect to volume: 


PvdW — — 


dTvdW 

dV 


= T 

- e 


’1-bp 




which at low density reduces to the ideal bulk swim pressure T^p. Going further, one can add a square gradient contri¬ 
bution to the free energy density, that simplifies the non local functional ( |34[ ) while allowing to describe inhomogeneous 
systems. This is achieved by introducing a term m|Vpp, where m = —{ljl2) ^ d^rr'^ We choose, now, 

the bare potential of the form w{r) = where Wq is the strength of the bare potential, cr is a nominal diameter 

and define f = rja, the non dimensional temperature T* = ^ and the parameter Kp = Pe^, where 

the last expression displays the dependence on the Peclet number. We, now, set : 


and 


Ts T* Pf.2a+2 J 


Ktjf) 


= - In 


l+2Ka 


1) 

, pOc -\-‘2 


1-2K, 


a 

Pra+l 


for f > 1. The effective hard-sphere diameter is given by the Barker-Henderson formula: 


CT Jo 


(36) 


(37) 


(38) 


and is density independent, but depends on the parameter Kp and temperature. As Kp —>■ 0 (i.e. at small Peclet 
number, i.e. small persistence time ) the system approaches the behavior of an assembly of passive soft spheres at 
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temperature T*, and inverse power law potentials analytic expressions for d{Ts) exist, whereas as Kp —> c» (large 
Peclet number) the activity becomes increasingly relevant and the effective radius increases with Kp. 

Let us remark that with the definition (371 of attractive potential the coefficient of formula (351 has a strong 


dependence on the temperature Tg, at variance with standard passive fluids where the dependence occurs through d{Tg) 
only. This makes the condensation transition basically athermal and driven by the activity parameter. This feature 
renders the case of active particles possessing also a direct attractive potential contribution of the Lennard-Jones 6-12 
type particularly interesting, that is obtained by modifying the bare pair potential as w\{r) = u>o[(f)“ — 

In fact, one has to consider a passive contribution similar to (351 to the vdW coefficient, say whose dependence 
on temperature is rather weak as compared to Thus one can expect that at low temperature and small Pe the 

condensation is driven by standard attraction, whereas at high Tg but large Pe the transition is determined by the 
effective attractive determined by the activitySl, 

We turn, now, to consider the virial series of the pressure: 


Tgp 


= {l + BiT)p + Cp^ + ...) 


and by comparing with the van der Waals equation we obtain second virial coefficient B{t) = {^pd^ — a^^\T)/Tg). 
Using the effective pair potential we can obtain B as the integral (with A = 0) 


B{T*,Kp) = 2Tra^ 

a{a + 1 ) 


1 + 2K, 


p 


pO;+2 


1 — exp 

rv n 2 


1 

T* \r 


1 


K 


Pp2a+2 I 


1 - 2K 


p fa+l 


}■ 


(39) 


Formula (391 gives the second virial coefficient for active spheres with purely repulsive potential a straightforward 


extension of it yields the values with finite values of A. 

The Boyle activity parameter, r* of the model, corresponds to the values where B = 0. For systems without 
attraction the curve Tg{T*) where B{Tg,T*) = 0 is monotonic and decreases when r* increases and the resulting 
phase diagram in a plane r, Tg is shown in Fig the region at the right of each line corresponds to i? < 0 for that 
given value of A. The same calculation is repeated for an active system with attraction and we found that the locus 
where B{Tg,T*) = 0 initially decreases as a function of r, but then it bends as displayed in Fig. and shows at 
sufficiently low temperatures Tg shows the presence of a region close to the origin where B < 0. Figure [^displays a 
reentrant behavior of the Boyle’s line for a = 12 when A > 0: at low temperatures T* the region B < 0 under the 
effect of the direct attractive force extends at small values of r, while the bending is totally absent in the case A = 0. 
The Barker-Henderson diameter diplays the same non monotonic trend along the Boyle line showing the correlation 
between the structural properties and the thermal properties of the system. 


VI. CONCLUSIONS 

In this paper we have investigated using a microscopic approach the steady state properties of a system of active 
particles and determined the stationary V-particle probability distribution function by requiring the vanishing of 
all currents. We used such a distribution to construct the BGY hierarchy for the reduce n-particle (with n < N) 
distribution functions and found that it contains terms which do not have a counterpart in systems of passive particles. 
The presence of these terms explains several experimental or numerical findings such the modified barometric law, the 
effective attraction between nominally repulsive objects, negative second virial coefficient, condensation. We defined 
an effective potential of attractive character due to the reduced mobility experienced by the active particles and 
determined by the combined effect of the excluded volume and the persistence of their motion. In an inhomogeneous 
environment the mobility turns out to be anisotropic and described by a tensor, whose expression we have obtained 
in the case of a spherical surface by separating its normal and tangential components. The same type of analysis 
allows to determine the mechanical pressure in systems of non interacting active particles via an hydrostatic balance, 
but the full extension to the interacting case remains an open problem. In order to obtain a bulk equation of state 
connecting pressure to density, swim temperature and activity parameter we have explored an alternative approach 
and constructed a mean field theory introducing a (pseudo) Helmholtz functional and determined the stationary 
distribution by a variational principle. In spite of the great similarities with equilibrium systems, profound differences 
remain. Already at the level of the first BGY equation one sees that it is not possible to disentangle the interaction 
terms stemming from the external potential from those from the particle-particle interactions. Such a feature in our 
opinion seems to obstruct the way to establish a full density functional theory. In fact, it does not seem possible to 
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xy 

Figure 4: Plane t,Ts'. the three curves in the main panel correspond to three different values (from right to left) of the 
attractive parameter A = 0,0.2, 0.4 and a = 12, respectively, and represent the loci where the second virial coefficient, B{t, Ts), 
vanishes. In the purely repulsive case (triangles) the behavior is monotonic, whereas when attraction is present the lines bend 
in the low temperature region, where B changes sign under the effect of the direct attractive force. In the inset we plot the 
effective Barker-Henderson diameter as a function of the effective temperature measured along the Boyle line in the three cases. 


prove the so called v-representability property, i.e. the unique correspondence between the one particle density profile 
and a given external potential, which is at the basis of DFT. On the hand, if one constructs the functional using 
the two body effective potential cj){r) the DFT approach sounds promising to tackle the properties of active particles 
under inhomogeneous condition. Finally we note that the homogeneous free energy that we have constructed seems 
to reproduce the behaviour of the active fluid, helps in locating its instabilities and the onset of phase separation. In 
this context it would be very interesting to compare our approximate theory with numerical simulations of interacting 
colored-noise driven particles. In the future we plan to apply the pseudo-Helmholtz functional to the analysis of 
strongly in-homogeneous systems of active particles and to investigate the MIPS in these situations. 
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Appendix A: Multidimensional unified colored noise approximation 


In the present appendix, following the method put forward by Cao and coworker d^^ l ^^ l, we introduce an auxiliary 
stochastic process, Wi, defined by: 


Wi = -h Vi. 

7 
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By taking the derivative of Wi with respect to t we obtain 


w. 


1 ■ 

= - } , T, - Xk + Vi- 

7 ^ dxk 
k 


(Al) 


after substituting into eq. (Al) relations (Qand (|^ to eliminate Vi and Xk we arrive at the evolution equation for Wi'. 




1 \" dFj I 1 

' k 


Wi - 

7 J 


D. 


1/2 


-’Htit)- 


(A2) 


Now, we assume the unified colored noise approximation (UCNAJP^ consisting in dropping the "acceleration term", 
Wi, featuring in the l.h.s. of this equation, so that we can write the following system of algebraic linear equations for 
the quantities Wi'. 


X ^ 

Oik -> — 

7 ^ OXk 


Wk 




7 V 


With the help of the matrix Tik (defined by eq. @) we go back to the equation for Xi eq.(Q , rewritten as Wi = 
Xi — to find the following Langevin equation: 

-* = E r.7 (f* + Dy^r.7,„(i) + D/yKi) + ^ ry ^ 

k ’ ' k I 

that we interpret using the Stratonovich conventiorP^. We finally observe that the last two terms can be gathered 
together and after some manipulations reported hereafter 



(A3) 

H eq.(0 , 

, rewritten as Wi = 

' dFk 
' dxi 

(A4) 


Eh = EEhl'r., + rr/I^],- = 


-iT dFk 




I k 


7 dxi 
■ 

7 dxi. 


kl 


7 dxi 




I k 


give rise to the form of the stochastic equation for Xi reported in eq. ©■ 

We limit ourselves to consider the stationary solutions with vanishing current and since we have assumed that 
the determinant is non vanishing the only solution of eq. (|^ corresponding to J; = 0 is obtained by imposing the 
conditions: 


1 


{Da + Dt)"f 


FkPN -PnYI = E 


dxi 


(A5) 


3 ■' 3 

Multiplying by T^k and summing over k after some algebra we obtain the first order differential equations for 


P: 


N 


{Da + Dt)^ 


E ^nkFk + Pn E gfc 

k jk 


dx. 


Bfen = 


dP, 


N 


dXr, 


Since the matrix elements contain 2-nd cross derivatives d'^jdxidxj of the potential U we have 


^ r - —r 

r kn — r\ ^ kj 
OXn 


dxj 


and we finally get 


d 


E g'fc ( ) “ E ^3k ( ] ~ 


jk \ J / j}^ 

The last equality derives from the following Jacobi’s formula 

1 d 


dXr, 


_ (j_ 

detr dxr, 


det r 


detr dy^. 


d 


(A6) 


(A7) 


(A8) 


(A9) 


where y stands for any of the variables Xi. 

Explicitly by using (A 61 and (A8) we find the system of differential equations ([^ determining the full probability 
distribution. 
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Appendix B: H-theorem for the FPE 


The functional + Ts Ih/at) discussed in Sect. IV has also an interesting dynamical property. In 

fact, after rewriting it under the form 

^:F[f^]=TrfN\n{^^-\nZN 


T. 


N 


where /nHx}, t) is a normalized solution of the dynamical Fokker-Planck equati on at i nstant t and P^r is the stationary 
distribution and the first term is the so called Kullback-Leibler relative entropjffiEZl^ One can show the following H- 
theorem 


dJ- 

dt 


< 0 . 


(Bl) 


The proof follows closely the method presented by RiskeiPSi (section 6.1 of his book), but before proceeding it is 
necessary to reduce the FPE described by eqs. ^ and (|^ to the canonical form: 


dfN 

dt 






(B2) 


where the diffusion matrix is 


and the drift vector is 


= (Da+ Dt)T-/ 


(B3) 


1 , 


.-1 

Now using some standard analytic manipulations reported by Risken it is simple to show the following formula 


^ + {Da + A)r-,^ ^F-, 


(B4) 


f = - / d-.f.D, 


.|( 2 ) 9Inii 9Ini? 
dxi dxj 


< 0 


(B5) 


where R = Jn/Pn- If is positive definite T must always decrease for yf 0 towards the minimum value 

—Ts IuZm- This result also implies that the solution of the FPE is unique, and after some time T the distance between 
two solutions is vanishingly small. 


Appendix C: Detailed Balance 


The detailed balance implies a stronger condition than the one represented by having a stationary distribution, 
since it implies that there is no net flow of probability around any closed cycle of states. In practice if detailed balance 
holds it is not possible to have a ratchet mechanism and directed motion. Again we use the equations derived by 
Risken (Section 6.4 of his bookwhich represent the sufficient and necessary conditions for detailed balan ce. Since 
the variables {xi} are even under time reversal we have to verify the validity of the following equationd^sBIl 


Pn{{x}) = Pn{{x}) 


dxj^ 


(Cl) 


Using the explicit form of (B31 and (B41 it is straightforward to verify that the detailed balance conditions are verified 


since they are the same as the condition expressing the vanishing of the current components Ji in the stationary state 
(see eq. (B2|. 


Appendix D: Approximation for the determinant 

The exact evaluation of the determinant F associated with the Hessian matrix is beyond the authors capabilities 
and we look for approximations in order to evaluate the effective forces. We consider the associated determinant in 
the case of two spatial dimensions: 














17 


/ [1 + 

-^w^^{r2,ri) 


T,j^i ^w^y{ri,rj) 

[1 + ^J2j^lWyy{ri,rj)] 

-^Wy^{r2,ri) 


-^'UJ^^{ri,r2) 

-^Wa^(ri,r2) 

[1+ ^Ej#2 ■w^=oir2,rj)] 


-^w^y{ri,rN) \ 

-^Wyy{ri,rN) 

-^w^y{r2,rN) 


\ -^w^y{rM,ri) 

and to order t/j is 


-^Wyy{rN,ri) 


[l + ^E,'^]vWya(riv,r,')] / 


detr 


1 + 


r 

7 


Y1 [wxxivi^rj) +Wyy{ri,rj)] 


(Dl) 


It is interesting to remark that the off-diagonal elements contain only one term, while the diagonal elements and their 
neighbors contain N elements. Thus in the limit of —)■ oo we expect that the matrix becomes effectively diagonal. 
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